### Setup R ###
###############

### Load packages
require(foreign)
require(sandwich)
require(lmtest)
require(MASS)
require(Hmisc)
require(stargazer)
require(sciplot)
require(lme4)

### Clear space
rm(list = ls())

### Load function to incorporate uncertainty (with robust standard errors)
milm <- function(fml, midata){
  xx <- terms(as.formula(fml))
  lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcovHC(tmp, type = "HC")))
    vcovs[[i]] <- vcovHC(tmp, type = "HC")
  }
  par.est <- apply(lms, 1, mean)
  se.within <- apply(ses, 1, mean)
  se.between <- apply(lms, 1, var)
  se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
  list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)   
}

### Setup Data ###
##################

### Load data. Data primarily from Fariss (2014), Keith et al (2009), and Keith (2012). See Appendix A for data source details and descriptions.
dat <- read.dta("supplemental.dta")

dat$fourfree <- as.numeric(dat$fourfree)
dat$PRESSLCK <- as.numeric(dat$PRESSLCK)
dat$STRIKELCK <- as.numeric(dat$STRIKELCK)
dat$HABEASLCK <- as.numeric(dat$HABEASLCK)
dat$PUBLICTRIALLCK <- as.numeric(dat$PUBLICTRIALLCK)
dat$FAIRTRIALLCK <- as.numeric(dat$FAIRTRIALLCK)
dat$TORTURELCK <- as.numeric(dat$TORTURELCK)
dat$TERMSLCK <- as.numeric(dat$TERMSLCK)
dat$FINALITYLCK <- as.numeric(dat$FINALITYLCK)
dat$EXCLUSIVEAUTHORITYLCK <- as.numeric(dat$EXCLUSIVEAUTHORITYLCK)
dat$BANEXCEPTLCK <- as.numeric(dat$BANEXCEPTLCK)
dat$FISCALAUTONOMYLCK <- as.numeric(dat$FISCALAUTONOMYLCK)
dat$SOPLCK <- as.numeric(dat$SOPLCK)
dat$ENUMQUALLCK <- as.numeric(dat$ENUMQUALLCK)
dat$JUDREVIEWLCK <- as.numeric(dat$JUDREVIEWLCK)
dat$HIERARCHLCK <- as.numeric(dat$HIERARCHLCK)
dat$LEGISLCK <- as.numeric(dat$LEGISLCK)
dat$DURATIONLCK <- as.numeric(dat$DURATIONLCK)
dat$DISSOLVELCK <- as.numeric(dat$DISSOLVELCK)
dat$LISTDEROLCK <- as.numeric(dat$LISTDEROLCK)
dat$COWCWAR2 <- as.numeric(dat$COWCWAR2)
dat$IW <- as.numeric(dat$IW)
dat$MIL2 <- as.numeric(dat$MIL2)
dat$LEFT <- as.numeric(dat$LEFT)
dat$BRIT <- as.numeric(dat$BRIT)
dat$POPTOTL  <- as.numeric(dat$POPTOTL)
dat$GDPPCPPPCURRENT <- as.numeric(dat$GDPPCPPPCURRENT)
dat$GLOBALNORMSICCPR <- as.numeric(dat$GLOBALNORMSICCPR)
dat$REGIONALNORMSICCPR  <- as.numeric(dat$REGIONALNORMSICCPR )

### Inspect data
colnames(dat)
attach(dat)

### Inspect values
summary(dat$PTSSDcorrect08)
summary(dat$l1PTSSDcorrect08)
summary(dat$PTSAIcorrected08)
summary(dat$l1PTSAIcorrected08)
summary(dat$fourfree)
summary(dat$PRESSLCK)
summary(dat$STRIKELCK)
summary(dat$HABEASLCK)
summary(dat$PUBLICTRIALLCK)
summary(dat$FAIRTRIALLCK)
summary(dat$TORTURELCK)
summary(dat$TERMSLCK)
summary(dat$FINALITYLCK)
summary(dat$EXCLUSIVEAUTHORITYLCK)
summary(dat$BANEXCEPTLCK)
summary(dat$FISCALAUTONOMYLCK)
summary(dat$SOPLCK)
summary(dat$ENUMQUALLCK)
summary(dat$JUDREVIEWLCK)
summary(dat$HIERARCHLCK)
summary(dat$LEGISLCK)
summary(dat$DURATIONLCK)
summary(dat$DISSOLVELCK)
summary(dat$LISTDEROLCK)
summary(dat$COWCWAR2)
summary(dat$IW)
summary(dat$MIL2)
summary(dat$LEFT)
summary(dat$BRIT)
summary(dat$latentmean)
summary(dat$LJI)
summary(dat$democracy)

### Create democracy dummy variables
dat$POLDEMO <- as.numeric(0)
dat$POLDEMO[DEMOC>5] <- as.numeric(1)
dat$POLDEMO[DEMOC==-88] <- NA
dat$POLDEMO[DEMOC==-77] <- NA
dat$POLDEMO[DEMOC==-66] <- NA
summary(dat$POLDEMO)
dat$FHDEMO <- as.numeric(0)
dat$FHDEMO[FHStatus>1] <- as.numeric(1)
dat$FHDEMO[FHStatus==-9] <- NA
summary(dat$FHDEMO)
dat$GWFDEM <- as.numeric(0)
dat$GWFDEM[gwf_nonautocracy=="democracy"] <- as.numeric(1)
summary(dat$GWFDEM)

attach(dat)

datfh <- na.omit(data.frame(latentmean,lag_latentmean,LJI,PRESSLCK, fourfree, STRIKELCK,HABEASLCK,PUBLICTRIALLCK,FAIRTRIALLCK,TORTURELCK,LEGISLCK,DURATIONLCK,DISSOLVELCK,LISTDEROLCK,COWCWAR2,IW,democracy,MIL2,LEFT,BRIT,POPTOTL,GDPPCPPPCURRENT,GLOBALNORMSICCPR,REGIONALNORMSICCPR,postsd,lag_latentsd,latentsd,uds_mean,uds_sd,COW,YEAR))

### Analyze Data ###
####################

### Estimate model with point estimates from latent variables
lmdd.out <- lm(latentmean ~ lag_latentmean + LJI + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR, data=datfh)
summary(lmdd.out)
lmdd.out$robustse <- vcovHC(lmdd.out)
a <- coeftest(lmdd.out,lmdd.out$robustse)
a

### Plot estimates
var.names <- c("De Facto Judicial Independence \n (Latent Measure)","Four Freedom Index","Freedom of Press","Right to Strike","Habeas Corpus","Public Trial","Fair Trial","Torture",
              "Legislative Declaration","Limited Duration","Can't Dissolve Legis.",
               "Non-Derogable Rights","Civil War","International War","Democracy","Military Control","State-Socialist Regime",
               "British Col. Exper.","Economic Development","Population","Global Norms","Regional Norms")
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.1, .1),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA)
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- coef(lmdd.out)[c(-1,-2)]
se <- sqrt(diag(vcovHC(lmdd.out)))[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#34495e")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#2980b9")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")

### Setup data to run models where we account for uncertainty
set.seed(12261991)
newdata <- list()
for(i in 1:1000){
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$latentmean, sd=datfh$latentsd)
  newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_latentmean, sd=datfh$lag_latentsd)
  newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJI, sd=datfh$postsd)
  newdata[[i]]$uds <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$uds_mean, sd=datfh$uds_sd)
}

### Check with uncertainty in DV and lagged DV
FML <- "draw ~ drawlag + LJI + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK +
FAIRTRIALLCK + TORTURELCK + LEGISLCK + DURATIONLCK + DISSOLVELCK + LISTDEROLCK + 
COWCWAR2 + IW + MIL2 + democracy + LEFT + BRIT + POPTOTL + GDPPCPPPCURRENT + 
GLOBALNORMSICCPR + REGIONALNORMSICCPR"
model <- milm(fml=FML, midata=newdata)
df <- 2231
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Plot estimates
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.75, .75),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- model$beta[c(-1,-2)]
se <- model$SE[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#2c3e50")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#d35400")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")

### Check with uncertainty in DV, lagged DV, and IV
FML2 <- "draw ~ drawlag + ljidraw + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + 
FAIRTRIALLCK + TORTURELCK + LEGISLCK + DURATIONLCK + DISSOLVELCK + LISTDEROLCK + 
COWCWAR2 + IW + MIL2 + democracy + LEFT + BRIT + POPTOTL + GDPPCPPPCURRENT + GLOBALNORMSICCPR + 
REGIONALNORMSICCPR"
model2 <- milm(fml=FML2, midata=newdata)
df <- 2231
tstats2<-data.frame(model2$terms,model2$beta,model2$SE,model2$beta/model2$SE,2*pt(-abs(model2$beta/model2$SE),df=df))
colnames(tstats2) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats2

### Plot estimates
plot(NULL,
     xlim = c(-.75, .75),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- model2$beta[c(-1,-2)]
se <- model2$SE[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#2c3e50")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#27ae60")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")

### Plot point estimates from all 3 models
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.5, .5),
     ylim = c(.9, 3),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 3) 
# without uncertainty
lines(c(0.1118128, 0.01555537), c(3, 3), lwd = 4, col="#95a5a6")
lines(c(0.103955, 0.02341311), c(3, 3), lwd = 10, col="#2980b9")
points(0.063684062285028, 3, pch = 19, cex = 1, col="white")
points(0.063684062285028, 3, pch = 19, cex = .5, col="black")
text(0.063684062285028, 3, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.063684062285028, 3, c("(Model 1)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in DV, and lagged DV
lines(c(0.5180008, 0.165589), c(2, 2), lwd = 4, col="#95a5a6")
lines(c(0.4892325, 0.1943573), c(2, 2), lwd = 10, col="#d35400")
points(0.3417949059079342, 2, pch = 19, cex = 1, col="white")
points(0.3417949059079342, 2, pch = 19, cex = .5, col="black")
text(0.3417949059079342, 2, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.3417949059079342, 2, c("(Model 2 - Uncertainty in DV, and Lagged DV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in DV, lagged DV, and IV
lines(c(0.4702866, 0.1414682), c(1, 1), lwd = 4, col="#95a5a6")
lines(c(0.4434443, 0.1683105), c(1, 1), lwd = 10, col="#27ae60")
points(0.3058773889594429, 1, pch = 19, cex = 1, col="white")
points(0.3058773889594429, 1, pch = 19, cex = .5, col="black")
text(0.3058773889594429, 1, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.3058773889594429, 1, c("(Model 3 - Uncertainty in DV, Lagged DV, and IV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 

### Cross Validation ###
########################

### Declare formulas and estimate OLS models
FORMULA <- "draw ~ drawlag"

model <- milm(fml=FORMULA, midata=newdata)
df <- 2252
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA2 <- "draw ~ drawlag + ljidraw"

model <- milm(fml=FORMULA2, midata=newdata)
df <- 2251
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA3 <- "draw ~ drawlag + democracy"

model <- milm(fml=FORMULA3, midata=newdata)
df <- 2251
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA4 <- "draw ~ drawlag + COWCWAR2"

model <- milm(fml=FORMULA4, midata=newdata)
df <- 2251
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA6 <- "draw ~ drawlag + fourfree + PRESSLCK + STRIKELCK + HABEASLCK +
PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LEGISLCK +
DURATIONLCK + DISSOLVELCK + LISTDEROLCK + COWCWAR2 +
IW + MIL2 + democracy + LEFT + BRIT + POPTOTL + GDPPCPPPCURRENT + 
GLOBALNORMSICCPR + REGIONALNORMSICCPR"

model <- milm(fml=FORMULA6, midata=newdata)
df <- 2232
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA7 <- "draw ~ drawlag + fourfree + PRESSLCK + STRIKELCK + HABEASLCK +
PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + ljidraw + LEGISLCK +
DURATIONLCK + DISSOLVELCK + LISTDEROLCK + COWCWAR2 +
IW + MIL2 + democracy + LEFT + BRIT + POPTOTL + GDPPCPPPCURRENT + 
GLOBALNORMSICCPR + REGIONALNORMSICCPR"

model <- milm(fml=FORMULA7, midata=newdata)
df <- 2231
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Declare functions
loss.rmse <- function(fit, df) cbind(as.numeric(predict(fit, newdata = df)), df$draw)

model.add <- function(df) lm(FORMULA, data=df)
model.add2 <- function(df) lm(FORMULA2, data=df)
model.add3 <- function(df) lm(FORMULA3, data=df)
model.add4 <- function(df) lm(FORMULA4, data=df)
model.add6 <- function(df) lm(FORMULA6, data=df)
model.add7 <- function(df) lm(FORMULA7, data=df)

validate.cv <- function(df, folds, resamples, model, loss) {
  
  lapply(1:resamples, function(x) {
    
    df$folds <- sample(rep(1:folds, length.out = nrow(df)))
    
    lapply(1:folds, function(test) {
      
      fit <- model(df[df$folds != test, ])
      
      loss(fit, df[df$folds == test, ])
      
    })})}

### Set values

SIMS <- 1000
cv1 <- NA
cv.add <- NA
cv2 <- NA
cv.add2 <- NA
cv3 <- NA
cv.add3 <- NA
cv4 <- NA
cv.add4 <- NA
cv5 <- NA
cv.add5 <- NA
cv6 <- NA
cv.add6 <- NA
cv7 <- NA
cv.add7 <- NA
newdata <- list()

for(i in 1:SIMS){
  
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$latentmean, sd=datfh$latentsd)
  newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_latentmean, sd=datfh$lag_latentsd)
  newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJI, sd=datfh$postsd)
  
  
  cv.add <- validate.cv(newdata[[i]], 10, 1, model.add, loss.rmse)
  cv.add2 <- validate.cv(newdata[[i]], 10, 1, model.add2, loss.rmse)
  cv.add3 <- validate.cv(newdata[[i]], 10, 1, model.add3, loss.rmse)
  cv.add4 <- validate.cv(newdata[[i]], 10, 1, model.add4, loss.rmse)  
  cv.add6 <- validate.cv(newdata[[i]], 10, 1, model.add6, loss.rmse)  
  cv.add7 <- validate.cv(newdata[[i]], 10, 1, model.add7, loss.rmse)  
  
  value <- rbind(cv.add[[1]][[1]],cv.add[[1]][[2]],cv.add[[1]][[3]],cv.add[[1]][[4]],cv.add[[1]][[5]],cv.add[[1]][[6]],cv.add[[1]][[7]],cv.add[[1]][[8]],cv.add[[1]][[9]],cv.add[[1]][[10]])
  cv1[i] <- sqrt(mean((value[,1] - value[,2])^2))
  value2 <- rbind(cv.add2[[1]][[1]],cv.add2[[1]][[2]],cv.add2[[1]][[3]],cv.add2[[1]][[4]],cv.add2[[1]][[5]],cv.add2[[1]][[6]],cv.add2[[1]][[7]],cv.add2[[1]][[8]],cv.add2[[1]][[9]],cv.add2[[1]][[10]])
  cv2[i] <- sqrt(mean((value2[,1] - value2[,2])^2))
  value3 <- rbind(cv.add3[[1]][[1]],cv.add3[[1]][[2]],cv.add3[[1]][[3]],cv.add3[[1]][[4]],cv.add3[[1]][[5]],cv.add3[[1]][[6]],cv.add3[[1]][[7]],cv.add3[[1]][[8]],cv.add3[[1]][[9]],cv.add3[[1]][[10]])
  cv3[i] <- sqrt(mean((value3[,1] - value3[,2])^2))
  value4 <- rbind(cv.add4[[1]][[1]],cv.add4[[1]][[2]],cv.add4[[1]][[3]],cv.add4[[1]][[4]],cv.add4[[1]][[5]],cv.add4[[1]][[6]],cv.add4[[1]][[7]],cv.add4[[1]][[8]],cv.add4[[1]][[9]],cv.add4[[1]][[10]])
  cv4[i] <- sqrt(mean((value4[,1] - value4[,2])^2))
  value6 <- rbind(cv.add6[[1]][[1]],cv.add6[[1]][[2]],cv.add6[[1]][[3]],cv.add6[[1]][[4]],cv.add6[[1]][[5]],cv.add6[[1]][[6]],cv.add6[[1]][[7]],cv.add6[[1]][[8]],cv.add6[[1]][[9]],cv.add6[[1]][[10]])
  cv6[i] <- sqrt(mean((value6[,1] - value6[,2])^2))
  value7 <- rbind(cv.add7[[1]][[1]],cv.add7[[1]][[2]],cv.add7[[1]][[3]],cv.add7[[1]][[4]],cv.add7[[1]][[5]],cv.add7[[1]][[6]],cv.add7[[1]][[7]],cv.add7[[1]][[8]],cv.add7[[1]][[9]],cv.add7[[1]][[10]])
  cv7[i] <- sqrt(mean((value7[,1] - value7[,2])^2))
}

y <- c(((cv2-cv1)/cv1),((cv3-cv1)/cv1),((cv4-cv1)/cv1),((cv6-cv1)/cv1),((cv7-cv1)/cv1))
a <- c(rep(c("(2) De Facto Judicial \n Independence \n Only"),1000),rep(c("(3) Democracy \n Only \n"),1000),rep(c("(4) Civil War \n Only \n"),1000),
       rep(c("(5) All Controls \n \n"),1000),rep(c("(6) Full Model \n \n "),1000))
d <- cbind(as.numeric(y),a)

### Plot reduction in mean square error across models
par(
  oma = c(0,0,0,0),
  mar = c(5,4,2,1)
)
bargraph.CI(d[,2],as.numeric(d[,1]), main="", col=c("#d35400","#cccccc","#cccccc","#cccccc","#d35400"), font.main=1,
            border=NA, family="Helvetica", ylim=c(min(as.numeric(d[,1])),0),err.width=.1,xlab="Models",ylab="Proportional Reduction in Mean Square Error")

### Robustness checks ###
#########################

### Appendix F, Model 2
set.seed(12261991)
newdata <- list()
for(i in 1:1000){
   newdata[[i]] <- datfh
   newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$latentmean, sd=datfh$latentsd)
   newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_latentmean, sd=datfh$lag_latentsd)
   newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJI, sd=datfh$postsd)
   newdata[[i]]$uds <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$uds_mean, sd=datfh$uds_sd)
 }

FML3 <- "draw ~ drawlag + ljidraw + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + 
FAIRTRIALLCK + TORTURELCK + LEGISLCK + DURATIONLCK + DISSOLVELCK + LISTDEROLCK + 
COWCWAR2 + IW + MIL2 + uds + LEFT + BRIT + POPTOTL + GDPPCPPPCURRENT + GLOBALNORMSICCPR + 
REGIONALNORMSICCPR"

model3 <- milm(fml=FML3, midata=newdata)
df <- 2231
tstats3<-data.frame(model3$terms,model3$beta,model3$SE,model3$beta/model3$SE,2*pt(-abs(model3$beta/model3$SE),df=df))
colnames(tstats3) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats3

### Unreported models

# base OLS model with Autocratic Regimes democracy measure
lmgwf.out <- lm(latentmean ~ lag_latentmean + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LJI + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + GWFDEM + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR, data=dat)
summary(lmgwf.out)
lmgwf.out$robustse <- vcovHC(lmgwf.out)
b <- coeftest(lmgwf.out,lmgwf.out$robustse)   

# base OLS model with Polity democracy measure
lmpol.out <- lm(latentmean ~ lag_latentmean + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LJI + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + POLDEMO + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR, data=dat)
summary(lmpol.out)
lmpol.out$robustse <- vcovHC(lmpol.out)
c <- coeftest(lmpol.out,lmpol.out$robustse)  

# base OLS model with Freedom House democracy measure
lmfh.out <- lm(latentmean ~ lag_latentmean + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LJI + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + FHDEMO + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR, data=dat)
summary(lmfh.out)
lmpol.out$robustse <- vcovHC(lmfh.out)
c <- coeftest(lmfh.out,lmfh.out$robustse)   

# base OLS model with State Department outcome measure
sdfh.out <- lm(PTSSDcorrect08 ~ l1PTSSDcorrect08 + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LJI + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR, data=dat)
summary(sdfh.out)
sdfh.out$robustse <- vcovHC(sdfh.out)
x <- coeftest(sdfh.out,sdfh.out$robustse)
   
# base OLS model with Amnesty International outcome measure
aifh.out <- lm(PTSAIcorrected08 ~ l1PTSAIcorrected08 + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LJI + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR, data=dat)
summary(aifh.out)
aifh.out$robustse <- vcovHC(aifh.out)
y <- coeftest(aifh.out,aifh.out$robustse)     

# country random effects
library(lme4)
lmrce.out <- lmer(latentmean ~ lag_latentmean + LJI + fourfree + PRESSLCK + STRIKELCK + HABEASLCK + PUBLICTRIALLCK + FAIRTRIALLCK + TORTURELCK + LEGISLCK + DURATIONLCK + 
                 DISSOLVELCK + LISTDEROLCK + COWCWAR2 + IW + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + POPTOTL + GLOBALNORMSICCPR + REGIONALNORMSICCPR + (1|COW), data=dat)
summary(lmrce.out)

# bivariate relationship
hrslji <- read.dta("hrslji.dta")
summary(lm(latentmean ~ LJI, data=hrslji))